In [1]:
from util import *
from glob import glob
/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:124: UserWarning: The Shapely GEOS version (3.11.2-CAPI-1.17.2) is incompatible with the GEOS version PyGEOS was compiled with (3.9.1-CAPI-1.14.2). Conversions between both will be slow.
  warnings.warn(
In [2]:
df = load_AOIs()
df
Out[2]:
Taranaki AOI SSP 4.5 (p50) SSP 4.5 (p83) SSP 8.5 (p50) SSP 8.5 (p83) Rate SSP 4.5 (p50) Rate SSP 4.5 (p83) Rate SSP 8.5 (p50) Rate SSP 8.5 (p83) match match_score
7 NORTH TongaporutuRiver 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 TongapurutuRiverCliffs 93.750000
11 SOUTH HangatahuaRiver_South 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 HangatahuRiver_South 97.435897
21 SOUTH Rahotu 0.58 0.78 0.84 1.10 0.0058 0.0078 0.0084 0.0110 Rahotu 100.000000
20 SOUTH Pihama 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Pihama 100.000000
19 SOUTH OpunakeBeach 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OpunakeBeachCliffs 100.000000
18 SOUTH OhaweBeach 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 OhaweBeach 100.000000
17 SOUTH Oeo 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Oeo 100.000000
16 SOUTH Manutahi 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Manutahi 100.000000
15 SOUTH ManaBay 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 ManaBayCliffs 100.000000
14 SOUTH KaupokonuiBeach 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 KaupokonuiBeach 100.000000
13 SOUTH Kakaramea 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Kakaramea 100.000000
12 SOUTH Hawera_WaihiBeach 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Hawera_WaihiBeach 100.000000
0 NORTH Mohakatino_River 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 MohakatinoRiver 100.000000
10 SOUTH CapeEgmont 0.58 0.78 0.84 1.10 0.0058 0.0078 0.0084 0.0110 CapeEgmont 100.000000
9 NORTH Waitara 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Waitara 100.000000
8 NORTH UrenuiRiver_North 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 UrenuiRiverNorthCliffs 100.000000
6 NORTH PariokariwaPoint 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 PariokariwaPointCliffs 100.000000
5 NORTH Onaero 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OnaeroCliff 100.000000
4 NORTH Oakura_South 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OakuraSouth 100.000000
3 NORTH Oakura 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Oakura 100.000000
2 NORTH NewPlymouth_Airport 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 NewPlymouthAirport 100.000000
1 NORTH NewPlymouth 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 NewPlymouthCliffs 100.000000
22 SOUTH WainuiBeach 0.57 0.78 0.83 1.09 0.0057 0.0078 0.0083 0.0109 WainuiBeach 100.000000
23 SOUTH WaipipiBeach 0.57 0.78 0.83 1.09 0.0057 0.0078 0.0083 0.0109 WaipipiBeachCliffs 100.000000
In [3]:
site = df.match.sample(1).iloc[0]
site
Out[3]:
'Manutahi'
In [4]:
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
if site == "ManaBayCliffs":
  print("Flipping")
  for k, v in transect_metadata.items():
    transect_metadata[k]["Azimuth"] = v["Azimuth"] + 180
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
In [5]:
results = predict(gdf, transect_metadata)
results
Out[5]:
TransectID BaselineID group Year ocean_point linear_model_point sqrt_model_point BH_model_point Sunamura_model_point
0 1 1 0 2100 POINT (1715590.6919486246 5609030.23268365) POINT (1715849.009481837 5609494.719162277) POINT (1715890.7038376718 5609569.6907071015) POINT (1715866.2782910739 5609525.770590698) POINT (1715854.3144944236 5609504.258222923)
1 2 1 0 2100 POINT (1715599.2111175938 5609023.965409) POINT (1715854.5722339319 5609484.838739857) POINT (1715894.3048528878 5609556.547795271) POINT (1715871.7922464146 5609515.917255742) POINT (1715859.8635021795 5609494.388370814)
2 3 1 0 2100 POINT (1715602.8857623506 5609014.427155222) POINT (1715859.3232794793 5609474.209777143) POINT (1715895.2309116768 5609538.590782948) POINT (1715876.6300122952 5609505.240084565) POINT (1715864.643914811 5609483.749472)
3 4 1 0 2100 POINT (1715601.6670005894 5609015.5665156795) POINT (1715868.964057981 5609471.345559046) POINT (1715908.6720163764 5609539.053207862) POINT (1715886.938245951 5609501.994075262) POINT (1715874.4882576733 5609480.76509578)
4 5 1 0 2100 POINT (1715601.965138068 5609016.699748106) POINT (1715879.6506328213 5609468.95492651) POINT (1715924.2571496856 5609541.603763499) POINT (1715898.2415991009 5609499.23328433) POINT (1715885.3620039735 5609478.256806616)
... ... ... ... ... ... ... ... ... ...
802 813 1 3 2100 POINT (1714827.4827839208 5609774.966452068) POINT (1715232.64006109 5610074.548259334) POINT (1715248.2960752654 5610086.1246453095) POINT (1715261.2087206184 5610095.672526711) POINT (1715241.452444554 5610081.064320906)
803 814 1 3 2100 POINT (1714833.7909443127 5609764.865698286) POINT (1715235.3934188318 5610064.101935908) POINT (1715243.5100513825 5610070.149683969) POINT (1715263.8844621596 5610085.330770769) POINT (1715244.186926914 5610070.654027693)
804 815 1 3 2100 POINT (1714841.9910656156 5609755.7839072) POINT (1715242.0664031254 5610056.6059552925) POINT (1715248.6121545136 5610061.527794142) POINT (1715270.4645472667 5610077.958903304) POINT (1715250.8322824126 5610063.197138285)
805 816 1 3 2100 POINT (1714851.789350054 5609747.469330703) POINT (1715250.379665398 5610050.344921249) POINT (1715258.1734205922 5610056.267137911) POINT (1715278.6693692633 5610071.841331058) POINT (1715259.1112112766 5610056.979734019)
806 817 1 3 2100 POINT (1714860.9832201214 5609738.160715329) POINT (1715257.424151746 5610043.13078496) POINT (1715265.2427249525 5610049.145377867) POINT (1715285.5857960042 5610064.794689783) POINT (1715266.1161317348 5610049.817263402)

807 rows × 9 columns

In [6]:
gpd.GeoSeries(results.linear_model_point).distance(gpd.GeoSeries(results.Sunamura_model_point)).describe()
Out[6]:
count    807.000000
mean      10.899348
std        0.022867
min       10.838591
25%       10.882344
50%       10.897024
75%       10.914473
max       10.967421
dtype: float64
In [7]:
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
poly = prediction_results_to_polygon(results)
In [8]:
poly.explore(tiles="Esri.WorldImagery")
Out[8]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [9]:
site = df.match.sample(1).iloc[0]
site = "CapeEgmont"
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
m = gpd.read_file(f"Shapefiles/{site}_TransectLines.shp").explore("TransectID", tiles="Esri.WorldImagery", max_zoom=22)
results = predict(gdf, transect_metadata)
results["geometry"] = results["linear_model_point"]
pd.concat([gdf, results])[["Year", "geometry"]].explore("Year", m=m)
gpd.GeoSeries(results.ocean_point, crs=2193).explore(m=m, color="blue")
#gpd.GeoSeries(poly, crs=2193).explore(color="pink", m=m)
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
polygon = prediction_results_to_polygon(results)
polygon.explore(m=m)
m
CapeEgmont
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
<ipython-input-9-91d19b3433a1>:10: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  results["geometry"] = results["linear_model_point"]
Out[9]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [10]:
run_all_parallel()
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
Flipping
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)

/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
In [11]:
def read_file(f):
  df = gpd.read_file(f)
  df["filename"] = f
  return df
samples = pd.concat(read_file(f) for f in glob("Projected_Shoreline_Polygons/*.shp"))
samples.explore("filename", tiles="Esri.WorldImagery", style_kwds={"fill": False})
Out[11]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [12]:
# Compare different SLR rates
site = df.match.sample(1).iloc[0]
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
low_proj_slr = predict(gdf, transect_metadata, Proj_SLR=0.007)
high_proj_slr = predict(gdf, transect_metadata, Proj_SLR=0.015)
m = gpd.GeoSeries(low_proj_slr.sqrt_model_point, crs=2193).explore(color="blue", tiles="Esri.WorldImagery")
gpd.GeoSeries(high_proj_slr.sqrt_model_point, crs=2193).explore(color="green", m=m)
Manutahi
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
Out[12]:
Make this Notebook Trusted to load map: File -> Trust Notebook